import iris
import matplotlib.pyplot as plt
/Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
Get some data...
import iris.fileformats.grib as grb
import iris.cube
import datetime
import cf_units
def dt_callback(cube, field, fname):
# Pull out the datetime from the filename.
loc_ECMF = fname.find('ECMF_')
loc_0h = fname.find('_0h')
dt_str = fname[loc_ECMF+5:loc_0h-2]
dt = datetime.datetime.strptime(dt_str, '%Y%m%d%H%M%S')
dt_unit = cf_units.Unit('days since 2000-01-01')
cube.add_aux_coord(iris.coords.DimCoord(standard_name='time',
points=[dt_unit.date2num(dt)],
units=dt_unit))
fname = 'A_HPXA89ECEM090000_C_ECMF_20160809000000_0h_em_msl_global_0p5deg_grib2.bin'
cubes = grb.load_cubes(fname, callback=dt_callback)
cubes = iris.cube.CubeList(cubes)
mslp = cubes.merge_cube()
print(mslp)
air_pressure / (Pa) (latitude: 361; longitude: 720) Dimension coordinates: latitude x - longitude - x Scalar coordinates: originating_centre: European Centre for Medium Range Weather Forecasts time: 2016-08-09 00:00:00 Attributes: GRIB_LOAD_WARNING: unsupported GRIB2 ProductDefinitionTemplate: #4.2
/Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/__init__.py:59: IrisDeprecation: The module iris.fileformats.grib is deprecated since v1.10. Please install the package 'iris_grib' package instead. "The module iris.fileformats.grib is deprecated since v1.10. " /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/__init__.py:935: IrisDeprecation: the`auto_regularise` kwarg is deprecated and will be removed in a future release. Resampling quasi-regular grids on load will no longer be available. Resampling should be done on the loaded cube instead using Cube.regrid. warn_deprecated(msg) /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/__init__.py:879: IrisDeprecation: Deprecated at version 1.10 warn_deprecated('Deprecated at version 1.10') /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/__init__.py:200: IrisDeprecation: Deprecated at version 1.10 warn_deprecated('Deprecated at version 1.10') /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/__init__.py:242: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future proxy = GribDataProxy(shape, np.zeros(.0).dtype, np.nan, /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/load_rules.py:69: IrisDeprecation: This GRIB loader is deprecated and will be removed in a future release. Please consider using the new GRIB loader by setting the :class:`iris.Future` option `strict_grib_load` to True; e.g.: iris.FUTURE.strict_grib_load = True Please report issues you experience to: https://groups.google.com/forum/#!topic/scitools-iris-dev/lMsOusKNfaU warn_deprecated(msg) /Users/pelson/miniconda/envs/iris_nb/lib/python2.7/site-packages/iris/fileformats/grib/load_rules.py:372: UserWarning: Different vertical bound types not yet handled. warnings.warn("Different vertical bound types not yet handled.")
Decimate the data (by a factor of 100), and then regrid the data onto this using a linear interpolation. Next, regrid the regridded data back up to the original grid, giving the overall effect of reducing the entropy.
mslp_smaller = mslp.regrid(mslp[::10, ::10], iris.analysis.Linear())
mslp = mslp_smaller.regrid(mslp, iris.analysis.Linear())
Get hold of the Australian coastline geometry
import cartopy.io.shapereader as shpreader
shpfilename = shpreader.natural_earth(resolution='50m',
category='cultural',
name='admin_0_countries')
reader = shpreader.Reader(shpfilename)
countries = reader.records()
[aus] = [country for country in countries if country.attributes['name'] == 'Australia']
aus.geometry
Produce a mask of points which lie within the geometry for our data's grid:
import shapely.vectorized
import numpy as np
x = mslp.coord('longitude').points
y = mslp.coord('latitude').points
x, y = np.meshgrid(x, y)
r = shapely.vectorized.contains(aus.geometry.buffer(0.15), x, y)
Mask the data outside of the geometry.
mslp.data[~r] = np.nan
Finally, plot the data...
import iris.plot as iplt
import cartopy.crs as ccrs
t = mslp.coord('time')
dt = t.units.num2date(t.points[0])
figsize = (7, 5)
with plt.xkcd():
plt.figure(figsize=figsize)
plt.axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
cs_f = iplt.contourf(mslp / 100, 6, corner_mask=True, cmap='viridis')
ax1 = plt.gca()
#cs = iplt.contour(mslp / 100, levels=cs_f.levels, colors=['black'], linewidths=[0.5])
plt.figure(figsize=figsize)
plt.axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
cs = iplt.contour(mslp / 100, levels=cs_f.levels, colors=['black'], linewidths=[0.5])
ax2 = plt.gca()
for ax in [ax1, ax2]:
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
#ax.coastlines('50m', zorder=10)
ax.add_geometries([aus.geometry.simplify(0.05, preserve_topology=True)], ax.projection, edgecolor='black', facecolor='none')
ax.set_extent([107, 165, -42, -7])
emblem_anchor_x, emblem_anchor_y = 116, -14.5
emblem = plt.imread('emblem.png')
ax.imshow(emblem, extent=[emblem_anchor_x - 4, emblem_anchor_x + 3,
emblem_anchor_y + 1.5, emblem_anchor_y + 6.5],
origin='upper')
ax.text(emblem_anchor_x, emblem_anchor_y, 'Australian Government\nBureau of Meteorology',
fontsize=10, ha='center', va='center', linespacing=2)
ax.plot([emblem_anchor_x - 7, emblem_anchor_x + 7],
[emblem_anchor_y, emblem_anchor_y],
color='black', lw=0.5)
ax.patch.set_visible(False)
ax.text(125, -40, 'Mean sea level pressure for {:%a %d, %b %Y}'.format(dt), fontsize=8, ha='center',
va='center', linespacing=2)
cbax = ax.figure.add_axes([0.85, 0.25, 0.035, 0.5])
cb = ax.figure.colorbar(cs_f, cax=cbax, label='hectopascal')
cb.ax.yaxis.label.set_fontsize(8)
ticks = cb.locator()
labels = np.arange(0,len(cs.levels),1)
cb.set_ticks(ticks[1:-1])
for i, tick in enumerate(ticks[:-1], start=1):
# Turn off the path effects for this piece of text so that we get color behind.
with plt.rc_context(rc={'path.effects': []}):
cb.ax.text(0.5, float(i - 0.5) / (len(ticks) - 1), i,
transform=cb.ax.transAxes, fontsize=7, va='center', ha='center')
cb.ax.tick_params(axis=u'both', which=u'both', length=0, labelsize=10)
from cartopy.mpl.patch import path_to_geos
for i, (level, collection) in enumerate(zip(cs_f.levels, cs_f.collections), start=1):
for path in collection._paths:
for geom in path_to_geos(path):
p = geom.representative_point()
if geom.area < 0.1:
# Don't bother with the little geometries (normally around the edge).
continue
for ax in [ax1, ax2]:
# Turn off the path effects for this piece of text so that we get color behind.
with plt.rc_context(rc={'path.effects': []}):
ax.text(p.x, p.y, i, fontsize=7, ha='center', va='center')
plt.show()